\chapter{Cost Function Framework}
\label{ch:cost_function}

The cost function framework in GSI represents the mathematical foundation for variational data assimilation. This chapter examines the implementation of the cost function evaluation and decomposition systems that enable optimal blending of observations with background forecasts through iterative minimization.

\section{Mathematical Foundation}

The GSI cost function represents a quadratic approximation to the maximum likelihood estimation problem in Bayesian statistics. The complete cost function is expressed as:

\begin{equation}
J(\mathbf{x}) = J_b(\mathbf{x}) + J_o(\mathbf{x}) + J_c(\mathbf{x}) + J_l(\mathbf{x}) + J_d(\mathbf{x})
\label{eq:total_cost_function}
\end{equation}

where:
\begin{itemize}
\item $J_b(\mathbf{x}) = \frac{1}{2} (\mathbf{x}-\mathbf{x}_b)^T \mathbf{B}^{-1} (\mathbf{x}-\mathbf{x}_b)$ is the background term
\item $J_o(\mathbf{x}) = \frac{1}{2} (\mathbf{y}-\mathbf{H}(\mathbf{x}))^T \mathbf{R}^{-1} (\mathbf{y}-\mathbf{H}(\mathbf{x}))$ is the observation term  
\item $J_c(\mathbf{x})$ represents constraint terms (dry mass, moisture)
\item $J_l(\mathbf{x})$ includes bias correction terms
\item $J_d(\mathbf{x})$ contains regularization terms
\end{itemize}

The control variable transformation enables preconditioning through:
\begin{equation}
\mathbf{x} = \mathbf{x}_b + \mathbf{B}^{1/2} \hat{\mathbf{x}}
\end{equation}

In hybrid ensemble-variational systems, the background error covariance becomes:
\begin{equation}
\mathbf{B} = (1-\alpha) \mathbf{B}_{\text{static}} + \alpha \mathbf{B}_{\text{ensemble}}
\end{equation}

\section{Cost Function Module Structure}

The \texttt{jfunc} module provides the core infrastructure for cost function minimization in GSI. Key components include:

\subsection{Control Vector Management}

\begin{lstlisting}[language=Fortran,caption={Cost Function Variables},label=code:jfunc_vars]
module jfunc
  ! Control variable dimensions
  integer(i_kind) :: nclen        ! Control vector length
  integer(i_kind) :: nval_len     ! Control subdomain size
  integer(i_kind) :: nvals_len    ! State subdomain size
  
  ! Iteration control
  integer(i_kind) :: jiter        ! Outer iteration counter
  integer(i_kind) :: jiterstart   ! First outer iteration
  integer(i_kind) :: jiterend     ! Last outer iteration
  integer(i_kind) :: niter        ! Inner iterations
  
  ! Solution storage
  type(control_vector) :: xhatsave ! Control variable storage
  type(control_vector) :: yhatsave ! Work vector storage
  
  ! Physical constraints
  logical :: clip_supersaturation
  real(r_kind) :: factqmin, factqmax
end module jfunc
\end{lstlisting}

The module manages memory allocation for large-scale control vectors with MPI-aware domain decomposition. Dynamic allocation scales with model resolution and observation types, supporting both global and regional configurations.

\subsection{Initialization System}

The \texttt{init\_jfunc} subroutine establishes default values for cost function parameters:

\begin{lstlisting}[language=Fortran,caption={Cost Function Initialization},label=code:init_jfunc]
subroutine init_jfunc()
  ! Outer iteration control
  first = .true.
  last = .false.
  jiterstart = 1
  jiterend = 1
  jiter = 1
  
  ! Minimization parameters  
  do i = 0, 30
    niter(i) = 0
    niter_no_qc(i) = 1000000
  end do
  
  ! Physical constraints
  factqmin = zero
  factqmax = zero
  clip_supersaturation = .false.
  
  ! Bias correction coefficients
  biascor(1) = 0.98_r_kind
  biascor(2) = 0.10_r_kind
  
  ! Algorithm controls
  switch_on_derivatives = .false.
  step_start = 1.e-4_r_kind
  iguess = 1
end subroutine init_jfunc
\end{lstlisting}

This initialization ensures consistent starting conditions while allowing flexibility through user-defined overrides in namelist configuration.

\section{Observation Cost Function Evaluation}

The \texttt{evaljo} subroutine computes the observation term of the cost function with high numerical precision:

\subsection{Mathematical Implementation}

The observation cost function is computed as:
\begin{equation}
J_o = \frac{1}{2} \sum_{i} w_i \frac{(y_i - H(x_b))^2}{R_i}
\end{equation}

where $w_i$ is the observation weight, $y_i$ is the observation value, $H(x_b)$ is the forward model prediction, and $R_i$ is the observation error variance.

\begin{lstlisting}[language=Fortran,caption={Observation Cost Function Evaluation},label=code:evaljo]
subroutine evaljo(pjo, kobs, kprt, louter)
  use kinds, only: r_quad
  use obsmod, only: nobs_type, cobstype, obscounts
  use mpl_allreducemod, only: mpl_allreduce
  
  real(r_quad), intent(out) :: pjo
  integer(i_kind), intent(out) :: kobs
  integer(i_kind), intent(in) :: kprt
  logical, intent(in) :: louter
  
  ! Local arrays for quadruple precision
  real(r_quad), dimension(nobs_bins*nobs_type) :: zprods
  real(r_quad), dimension(nobs_bins, nobs_type) :: zjo2
  real(r_quad), dimension(nobs_type) :: zjo1
  
  ! Initialize contributions
  zprods = zero_quad
  iobsgrp = 0
  
  ! Loop over observation types and bins
  jj = 0
  do jtype = 1, nobs_type
    do jbin = 1, nobs_bins
      jj = jj + 1
      
      ! Process observations in this type/bin
      obsptr => obsdiags(jtype, jbin)%ptr
      do while (associated(obsptr))
        if (obsptr%luse .and. obsptr%muse(jiter)) then
          ! Compute departure
          if (louter) then
            zdep = obsptr%nldepart
          else
            zdep = obsptr%tldepart - obsptr%nldepart  
          end if
          
          ! Accumulate weighted square
          zprods(jj) = zprods(jj) + obsptr%wgtjo * zdep**2
        end if
        obsptr => obsptr%next
      end do
    end do
  end do
  
  ! Global reduction
  call mpl_allreduce(nobs_bins*nobs_type, zprods)
  call mpi_allreduce(iobsgrp, iobsglb, nobs_bins*nobs_type, &
                    mpi_integer, mpi_sum, mpi_comm_world, ierror)
  
  ! Final aggregation
  pjo = zero_quad
  kobs = 0
  do jtype = 1, nobs_type
    do jbin = 1, nobs_bins
      jj = (jtype-1)*nobs_bins + jbin
      zjo2(jbin,jtype) = zprods(jj)
      pjo = pjo + zjo2(jbin,jtype)
      kobs = kobs + iobsglb(jbin,jtype)
    end do
  end do
end subroutine evaljo
\end{lstlisting}

\subsection{Precision and Reproducibility}

The use of quadruple precision arithmetic ensures:

\begin{itemize}
\item Reproducible results across different computer architectures
\item Numerical stability during iterative minimization
\item Accurate aggregation of contributions from millions of observations
\item Consistent behavior in parallel processing environments
\end{itemize}

\section{Cost Function Decomposition}

The GSI system provides detailed decomposition of cost function components for diagnostic purposes:

\subsection{Per-Observation Type Analysis}

\begin{lstlisting}[language=Fortran,caption={Cost Function Diagnostics},label=code:jo_diagnostics]
if (kprt >= 2) then
  write(iout_jo,*) 'Cost function Jo analysis:'
  write(iout_jo,'(a10,20a12)') 'Type\\Bin', (jbin,jbin=1,nobs_bins)
  
  do jtype = 1, nobs_type
    write(iout_jo,'(a10,20f12.4)') trim(cobstype(jtype)), &
          (real(zjo2(jbin,jtype)/max(iobsglb(jbin,jtype),1)), jbin=1,nobs_bins)
  end do
  
  write(iout_jo,*) 
  write(iout_jo,'(a,f15.6)') 'Total Jo:', real(pjo)
  write(iout_jo,'(a,i8)') 'Total observations:', kobs
  write(iout_jo,'(a,f15.6)') 'Jo per observation:', real(pjo)/max(kobs,1)
end if
\end{lstlisting}

This diagnostic output enables:
\begin{itemize}
\item Identification of observation types with poor fit
\item Detection of systematic biases in specific observation categories
\item Monitoring convergence behavior by observation type
\item Quality control validation through cost function analysis
\end{itemize}

\section{Quality Limit Evaluation}

The \texttt{evalqlim} subroutine implements quality control limits for specific observation types:

\begin{lstlisting}[language=Fortran,caption={Quality Limit Implementation},label=code:evalqlim]
subroutine evalqlim(ges_q, ges_qsat, qlim_adj)
  ! Apply moisture quality control limits
  ! Prevent supersaturation in analysis
  
  do k = 1, nsig
    do j = 1, lon2
      do i = 1, lat2
        ! Compute saturation mixing ratio
        call genqsat(qsat, ges_tsen(i,j,k), ges_prsl(i,j,k))
        
        ! Apply supersaturation constraint
        if (clip_supersaturation) then
          qlim = factqmax * qsat
          if (ges_q(i,j,k) > qlim) then
            qlim_adj(i,j,k) = qlim - ges_q(i,j,k)
          end if
        end if
        
        ! Apply minimum humidity constraint  
        if (ges_q(i,j,k) < factqmin * qsat) then
          qlim_adj(i,j,k) = factqmin * qsat - ges_q(i,j,k)
        end if
      end do
    end do
  end do
end subroutine evalqlim
\end{lstlisting}

\section{Hybrid Ensemble Integration}

The cost function framework seamlessly integrates with hybrid ensemble-variational methods:

\subsection{Ensemble Background Error}

\begin{lstlisting}[language=Fortran,caption={Hybrid Ensemble Cost Function},label=code:hybrid_cost]
if (l_hyb_ens) then
  ! Hybrid ensemble-variational cost function
  ! B = (1-alpha) * B_static + alpha * B_ensemble
  
  ! Static component
  call evaljo_static(pjo_static, kobs_static)
  
  ! Ensemble component  
  call evaljo_ensemble(pjo_ensemble, kobs_ensemble)
  
  ! Weighted combination
  pjo = (1.0_r_kind - ensemble_alpha) * pjo_static + &
        ensemble_alpha * pjo_ensemble
        
  ! Localization application
  if (l_localization) then
    call apply_localization(pjo_ensemble)
  end if
end if
\end{lstlisting}

\subsection{Dual Resolution Support}

For computational efficiency, the ensemble component may be computed on a coarser grid:

\begin{lstlisting}[language=Fortran,caption={Dual Resolution Implementation},label=code:dual_res]
if (dual_res) then
  ! Interpolate ensemble increments to analysis grid
  call ensemble_forward_model(xhat_ens, sval_ens)
  call interpolate_ens2anl(sval_ens, sval_interp)
  
  ! Apply ensemble background error on analysis grid
  call evaljo_dual(sval_interp, pjo_ens)
else
  ! Single resolution - direct evaluation
  call evaljo_single(xhat, pjo)
end if
\end{lstlisting}

\section{Memory Management and Performance}

The cost function framework implements several optimization strategies:

\subsection{Memory Allocation Strategy}

\begin{lstlisting}[language=Fortran,caption={Dynamic Memory Management},label=code:memory_mgmt]
subroutine create_jfunc()
  ! Allocate control vectors based on configuration
  call allocate_cv(xhatsave)
  call allocate_cv(yhatsave)
  
  ! State vector allocation
  call allocate_state(sval)
  
  ! Bias predictor allocation
  if (nsclen > 0) call allocate_preds(bval)
  
  ! Observation diagnostics
  allocate(obscounts(nobs_bins, nobs_type))
end subroutine create_jfunc

subroutine destroy_jfunc()
  ! Clean deallocation
  call deallocate_cv(xhatsave)
  call deallocate_cv(yhatsave)
  call deallocate_state(sval)
  if (allocated(obscounts)) deallocate(obscounts)
end subroutine destroy_jfunc
\end{lstlisting}

\subsection{Parallel Processing Considerations}

The cost function evaluation leverages MPI parallelization:

\begin{itemize}
\item \textbf{Domain decomposition}: Observations distributed across processors
\item \textbf{Collective operations}: Global reductions for cost function summation
\item \textbf{Load balancing}: Dynamic work distribution based on observation density
\item \textbf{Communication optimization}: Minimized data movement through local computation
\end{itemize}

\section{Diagnostic Capabilities}

The cost function framework provides comprehensive monitoring capabilities:

\subsection{Convergence Monitoring}

\begin{lstlisting}[language=Fortran,caption={Convergence Diagnostics},label=code:convergence]
! Track cost function reduction
if (jiter > 1) then
  cost_reduction = (pjo_prev - pjo) / pjo_prev
  write(iout_iter,'(a,i3,2f15.6,f10.4)') &
    'Iter:', jiter, pjo_prev, pjo, cost_reduction
  
  ! Check convergence criteria
  if (cost_reduction < convergence_threshold) then
    converged = .true.
  end if
end if
\end{lstlisting}

\subsection{Quality Control Integration}

The cost function framework integrates with the GSI quality control system:

\begin{itemize}
\item \textbf{Nonlinear QC}: Iterative quality control based on cost function analysis
\item \textbf{Variational QC}: Simultaneous quality control and analysis
\item \textbf{Adaptive observation weighting}: Dynamic weight adjustment based on fit
\item \textbf{Outlier detection}: Statistical identification of problematic observations
\end{itemize}

\section{Summary}

The cost function framework in GSI provides a robust and flexible foundation for variational data assimilation. Key features include:

\begin{itemize}
\item Mathematical rigor through proper quadratic cost function formulation
\item High numerical precision ensuring reproducible results
\item Comprehensive diagnostic capabilities for monitoring and validation
\item Seamless integration with hybrid ensemble-variational methods
\item Efficient parallel implementation with optimized memory management
\item Flexible quality control and constraint application
\end{itemize}

This framework enables GSI to achieve optimal analysis quality while maintaining computational efficiency and numerical stability across diverse atmospheric prediction applications.